Fast algorithm for calculating two-photon absorption spectra 
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We report a numerical calculation of the two-photon absorption coefficient of electrons in a binding 
potential using the real-time real-space higher-order difference method. By introducing random 
vector averaging for the intermediate state, the task of evaluating the two-dimensional time integral 
is reduced to calculating two one-dimensional integrals. This allows the reduction of the computation 
load down to the same order as that for the linear response function. The relative advantage of 
the method compared to the straightforward multi-dimensional time integration is greater for the 
calculation of non-linear response functions of higher order at higher energy resolution. 
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The measurement of two-photon absorption coefficient 
yields different information from the single-photon ab- 
sorption measurement, since the physical processes in- 
volved and selection rules are different. Because of this, 
there has always been a lot of interest in two-photon ab- 
sorption of various molecules, crystals and solids 
To compare against the experimental data, one would 
like to have theorerical calculations based on some re- 
alistic model of the material. However, since the two- 
photon absorption is a typical non-linear optical process, 
its realistic modelinghas always proved difficult for large 
complex systems []2j-|6|. 

A powerful method that has come to be used widely for 
large quantum systems is the real-time real-space higher- 
order difference method ]X0| , |Tl[] , in which the real-space 
is represented by discrete mesh points, and the time 
development of a system is solved by numerically inte- 
grating the Schrodinger equation for discrete time steps. 
The energy levels and energy eigenstates are obtained 
by Fourier analyzing the numerical solution. The mem- 
ory requirement scales linearly with the number of basis 
states N compared to iV 2 for matrix diagonalization, and 
the method has proved effective in solving large quantum 
systems that cannot be solved by conventional methods 
]ll],[t2| . So far, the large computation load has meant 
that the actual application of the method has been made 
primarily to the calculation of linear response functions 
of one-particle systems |lC],[ll|. Nevertheless, the poten- 
tial scale advantage of the method when applied to large 
systems invites the speculation that the development of 
the method will be essential in making the calculation 
of non-linear response functions of complex many-body 
quantum systems feasible. 

In this article, we report a trial application of the 
method to the calculation of the two-photon absorption 
coefficient of non-interacting electrons trapped by a bind- 
ing potential and exposed to monochromatic light of fre- 
quency u>. 

The two-photon absorption coefficient is given by |IJ 
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where H is the unperturbed hamiltonian of the system, 
r is the electron coordinate operator, r) is the frequency 
resolution, V is the volume of the system, the summation 
with respect to the initial state \Ei) must be taken over 
all states below the Fermi level E Fl and the summation 
with respect to \Ef) over the states Ef > E F - 

An important ingredient of the real-time real-space 
higher-order difference method is the use of random vec- 
tors as probes to scan the Hilbert space ]lO| , |TT| , |T3| [L6| . 
Among various types of random vector is the uniform 
amplitude random phase vector 
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whose effectiveness has been amply demonstrated 
|pT^ , p]Jl3| -[rrj] . Here, the phases 4> n are independent ran- 
dom variables with uniform distribution in the range 
[— 7r, 7r), and \n) (n = 1,N) are the orthonormal basis 
states which are localized at the mesh points in the real 
space. Using the property of completeness 
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which obtains after averaging over random realizations of 
|<E>) denoted above by the brackets (• • ■)«, the two-photon 
absorption coefficient eq.(fi]) can be rewritten as 
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The operator step function j[5) 
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can be explicitly constructed for any bounded hermi- 
tian operator X without solving for the eigenvalues Xi 
and eigenvectors \Xi). In our calculation, we used an 
algorithm based on Chebyshev polynomial expansion, 
which yields 6{X) as a polynomial of the operator X 

If eq.(|4|) were to be numerically implemented straight- 
forwardly, the matrix element inside the integral would 
have to be obtained for all necessary combinations of time 
variables t\ and ■ One would then start with a random 
vector |<fr), solve the time development according to the 
Schrodinger equation, mutiplying operators (r's and 0's) 
along the way as required, and take the inner product 
with |$) at the end. The number of discrete time steps 
required for a calculation of energy resolution ri scales as 
77 -1 , so that the direct implementation of eq.(0) requires 
a computation load that grows as rj~ 2 . On top of it, 
the calculation has to be repeated for a number of dif- 
ferent realizations of |3>) for random averaging. This is 
necessary to reduce the fluctuation arising from the use 
of random vectors, whose amplitude can be of the same 
order in magnitude as the final result itself. The scale of 
such computation can easily overwhelm the capacity of 
any computing facility in existence. 

However, the computational load can be greatly re- 
duced by inserting the completeness relation eq.(|3j) in 
the matrix element of eq.(^) to decompose it into two 
factors. The absorption coefficient is then given by 
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where |$) and |$') are mutually independent random 
vectors, and E c is the cutoff energy to be explained below. 
The most costly process of integrating the Schrodinger 
equation is now used only to obtain two complex func- 
tions of a single time variable instead of a bivariate 
function with two time variables. Once the neces- 
sary complex-valued functions have been calculated and 
stored, the two-dimensional time integration may easily 
be done with a small computer. The computational load 
therefore scales only as r\~ x with the energy resolution. 

The benefit must be weighed against the increased cost 
of having to average over random realizations of interme- 
diate states |$'). The statistical variance arising from the 
random sampling is independent of the number of time 
steps used in the calculation, but is controlled only by 
the number of random samples taken. Therefore, the 
relative advantage of the use of eq. @ over the straight- 
forward integration of eq.(Q) increases as higher energy 
resolution is required. For the actual numerical imple- 
mentation of the random sampling, it is essential that 



one has control over the extent of the Hilbert space to 
be probed. In eq.(^), the extra cutoff factor 9(E C — H) 
is inserted for this purpose. The final result should be 
independent of the cutoff energy E c if it is taken suffi- 
ciently large. However, a large value of E c entails large 
random fluctuation, so that larger number of samples will 
be required to average out the noise to achieve the same 
accuracy. Therefore, the value of E c must be set as small 
as possible provided that its interference with the final 
result is kept within the margin of tolerance. 

In order to compare the CPU time for the two meth- 
ods, we have computed eqs.(||) and (||) for the case of a 
parabolic binding potential 
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with wo = 0.3 atomic unit (a.u.) and m being the 
electron mass. The real space was represented by 16 3 
mesh points to cover a cubic volume of linear dimen- 
sion 16 a.u. For the time development, discrete time 
steps with At = 0.05 a.u. were used for integration of 
the Schrodinger equation. The Fermi energy was set at 
Ep = 3Tiu!o and the frequency resolution was rj = 8 x 10~ 2 
a.u. 
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FIG. 1. The two-photon absorption coefficient (uj) 
of electrons trapped in a parabolic potential. The results 
obtained by the use of eqs.(4) and (6) are shown by the dotted 
curve and the solid curve, respectively. The analytical result 
is shown by the dashed line. 

In fig.Qwe compare the results of eqs.(Q) and (||) with 
the analytical result for electrons in the parabolic poten- 
tial. Both numerical results well reproduce the analytical 
curve with the standard deviation standing at 7% at the 
peak (u>o = 0.3 a.u.) for both cases. 

The result for the straightforward implementation of 
eq.fl) is an average over 200 runs with different |$)'s. 
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The CPU time on a single processing unit of Fujitsu 
VPP500 was 1.5 x 10 3 seconds for each run, totaling 
3.0 x 10 5 seconds (~ 83 hours) to achieve the 7% accu- 
racy. For the calculation according to eq.(|J), an average 
was first taken over 50 different samples of |$') with a 
fixed |$). The result was then averaged over 100 different 
samples of |$). The total of 5000 runs of integrating the 
Schrodinger equation took 6.8 x 10 4 seconds (~ 19 hours), 
to achieve the same 7% accuracy at the peak. If the sta- 
tistical standard deviation is to be brought down to 1%, 
an average over 7 2 = 49 times more samples will have to 
be taken, which translates to 1.5 x 10 7 seconds (~ 174 
days) and 3.3 x 10 6 seconds (~ 38 days) for eqs.(^) and 
(0), respectively. In fig. 2, we show the estimated CPU 
time to achieve 1% statistical accuracy for various values 
of frequency resolution. The relative advantage of the 
use of random vectors for the intermediate state should 
grow as higher frequency resolution is required. The ac- 
tual lapse-time of computation can be reduced nearly by 
an order of magnitude if the computation is parallelized 
to use all the 30 processing units on Fujitsu VPP500 at 
RIKEN. 
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FIG. 2. The relation between the energy resolution and 
the total CPU time for calculating the two-photon absorption 
coefficient to 1% statistical accuracy. The solid line and the 
dotted line are for calculation according to eqs.(6) and (4) of 
the text, respectively. The cross and the circle are the pro- 
jected CPU time for energy resolution r\ — 8 x 1CP 2 a.u. based 
on the actual calculation performed with the same resolution 
but to 7% accuracy. 



tion to electrons in various pseudo-potentials. Extension 
of the algorithm itself to deal with an interacting M- 
particle system is straightforward, although solution of 
the Schrodinger equation in the 3M-dimensional space 
soon becomes unmanageable with increasing M. Of more 
significance is the potential of the present strategy to be 
used to facilitate the calculation of nonlinear response 
functions of higher order. A numerical calculation of a 
non-linear response function of order n typically involves 
computing an n-dimensional integral of the type 
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where each O represents some operator. While the 
straightforward evaluation requires computation time 
proportional to r/ _Tl , it may be possible to reduce the 
CPU time as far down as to n times that for linear re- 
sponse function by inserting the completeness relation 
cq.(||) and decomposing the matrix element into n fac- 
tors. This may be regarded as a version of quantum 
Monte Carlo method, and the use of importance sam- 
pling techniques (l4|jll| will be vital in reducing the eval- 
uation time. In fact, the use of 9(E C — H) in eq.(||) is 
one rudimentary method of improving the sampling effi- 
ciency. 

The insertion of random intermediate vectors is not 
the only way to reduce the CPU time of evaluation 
of the n-dimcnsional integral. For example, the stan- 
dard Monte Carlo random sampling technique may 
be applied directly to the evaluation of the integral 
eq.(0). If one adopts uniform random sampling in the 
n-dimensional time space, the CPU time for integration 
of the Schrodinger equation scales as r/ _1 times the num- 
ber of samples required for statistical averaging, which 
is precisely the same as for the case of random vector 
insertion above. In cither case, the statistical method of 
evaluation suffers from the familiar negative sign prob- 
lem jl!|, so that some scheme needs to be devised to 
improve the sampling efficiency. Nevertheless, it remains 
true that the size of statistical variance is independent 
of the energy resolution, so that it should be advanta- 
geous to employ statistical methods for the calculation 
of high-order nonlinear coefficient at high resolution of 
frequency. 

All calculations reported here were performed on a Fu- 
jitsu VPP500 at RIKEN. 



The present calculation was done for the specific case 
of non-interacting electrons trapped in a parabolic po- 
tential, only in order to test the relative advantage of 
the algorithm and to compare against the analytical re- 
sult. The computer code admits an arbitrary potential, 
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